%%%% A separate Matlab code to plot transition dynamics

%% initilization
start_time = 0; end_time = 30;
time_period = start_time: (end_time + 1);

K_low_baseline     = zeros(length(time_period),1); K_high_baseline = K_low_baseline; K_low_eta = K_low_baseline; K_high_eta = K_low_baseline;
B_low_baseline     = zeros(length(time_period),1); B_high_baseline = B_low_baseline; B_low_eta = B_low_baseline; B_high_eta = B_low_baseline;
Y_low_baseline     = zeros(length(time_period),1); Y_high_baseline = B_low_baseline; Y_low_eta = B_low_baseline; Y_high_eta = B_low_baseline;
I_low_baseline     = zeros(length(time_period),1); I_high_baseline = B_low_baseline; I_low_eta = B_low_baseline; I_high_eta = B_low_baseline;
B_Y_low_baseline   = zeros(length(time_period),1); B_Y_high_baseline = B_low_baseline; B_Y_low_eta = B_low_baseline; B_Y_high_eta = B_low_baseline;
C_low_baseline     = zeros(length(time_period),1); C_high_baseline = B_low_baseline; C_low_eta = B_low_baseline; C_high_eta = B_low_baseline;
L_low_baseline     = zeros(length(time_period),1); L_high_baseline = B_low_baseline; L_low_eta = B_low_baseline; L_high_eta = B_low_baseline;
phi_low_baseline   = zeros(length(time_period),1); phi_high_baseline = B_low_baseline; phi_low_eta = B_low_baseline; phi_high_eta = B_low_baseline;
q_low_baseline     = zeros(length(time_period),1); q_high_baseline = B_low_baseline; q_low_eta = B_low_baseline; q_high_eta = B_low_baseline;
tau_l_low_baseline = zeros(length(time_period),1); tau_l_high_baseline = B_low_baseline; tau_l_low_eta = B_low_baseline; tau_l_high_eta = B_low_baseline;
tau_k_low_baseline = zeros(length(time_period),1); tau_k_high_baseline = B_low_baseline; tau_k_low_eta = B_low_baseline; tau_k_high_eta = B_low_baseline;
R_low_baseline     = zeros(length(time_period),1); R_high_baseline = B_low_baseline; R_low_eta = B_low_baseline; R_high_eta = B_low_baseline;
Gov_spending_low_baseline     = zeros(length(time_period),1); Gov_spending_high_baseline = B_low_baseline; Gov_spending_low_eta = B_low_baseline; Gov_spending_high_eta = B_low_baseline;

q_w_low_baseline     = zeros(length(time_period),1); q_w_high_baseline = B_low_baseline; q_w_low_eta = B_low_baseline; q_w_high_eta = B_low_baseline;
r_K_low_baseline     = zeros(length(time_period),1); r_K_high_baseline = B_low_baseline; r_K_low_eta = B_low_baseline; r_K_high_eta = B_low_baseline;
rho_low_baseline     = zeros(length(time_period),1); rho_high_baseline = B_low_baseline; rho_low_eta = B_low_baseline; rho_high_eta = B_low_baseline;

%% Load AMPL results and transform the results to Matlab

%%%% transition, traditional Ramsey model
load transition_baseline_eta0
for t = start_time:end_time+1 
    K_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.K',num2str(t)));
    B_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.B',num2str(t)));    
    Y_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.Y',num2str(t)));    
    B_Y_low_baseline(t-start_time+1)   = eval(strcat('transition_traditional.B_Y',num2str(t)));
    C_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.C',num2str(t)));
    L_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.L',num2str(t)));
    phi_low_baseline(t-start_time+1)   = eval(strcat('transition_traditional.phi',num2str(t)));
    I_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.I',num2str(t)));
    q_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.q',num2str(t)));
    tau_l_low_baseline(t-start_time+1) = eval(strcat('transition_traditional.tau_l',num2str(t)));
    tau_k_low_baseline(t-start_time+1) = eval(strcat('transition_traditional.tau_k',num2str(t)));
    R_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.R',num2str(t)));
    Gov_spending_low_baseline(t-start_time+1)     = eval(strcat('transition_traditional.G',num2str(t)));  
end

%%%% transition, Ramsey model with financial frictions (our main model)
load transition_baseline_eta_positive
for t = start_time:end_time+1 
    K_high_baseline(t-start_time+1)     = eval(strcat('transition.K',num2str(t)));
    B_high_baseline(t-start_time+1)     = eval(strcat('transition.B',num2str(t)));    
    Y_high_baseline(t-start_time+1)     = eval(strcat('transition.Y',num2str(t)));    
    B_Y_high_baseline(t-start_time+1)   = eval(strcat('transition.B_Y',num2str(t)));
    C_high_baseline(t-start_time+1)     = eval(strcat('transition.C',num2str(t)));
    L_high_baseline(t-start_time+1)     = eval(strcat('transition.L',num2str(t)));
    phi_high_baseline(t-start_time+1)   = eval(strcat('transition.phi',num2str(t)));
    I_high_baseline(t-start_time+1)     = eval(strcat('transition.I',num2str(t)));
    q_high_baseline(t-start_time+1)     = eval(strcat('transition.q',num2str(t)));
    tau_l_high_baseline(t-start_time+1) = eval(strcat('transition.tau_l',num2str(t)));
    tau_k_high_baseline(t-start_time+1) = eval(strcat('transition.tau_k',num2str(t)));
    R_high_baseline(t-start_time+1)     = eval(strcat('transition.R',num2str(t)));
    Gov_spending_high_baseline(t-start_time+1)     = eval(strcat('transition.G',num2str(t))); 
end

%%%% financial shocks, traditional Ramsey model
load transition_fiscal_shocks_eta0
for t = start_time:end_time+1 
    K_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.K',num2str(t)));
    B_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.B',num2str(t)));    
    Y_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.Y',num2str(t)));    
    B_Y_low_eta(t-start_time+1)   = eval(strcat('transition_traditional_fiscal_shocks.B_Y',num2str(t)));
    C_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.C',num2str(t)));
    L_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.L',num2str(t)));
    phi_low_eta(t-start_time+1)   = eval(strcat('transition_traditional_fiscal_shocks.phi',num2str(t)));
    I_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.I',num2str(t)));
    q_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.q',num2str(t)));
    tau_l_low_eta(t-start_time+1) = eval(strcat('transition_traditional_fiscal_shocks.tau_l',num2str(t)));
    tau_k_low_eta(t-start_time+1) = eval(strcat('transition_traditional_fiscal_shocks.tau_k',num2str(t)));
    R_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.R',num2str(t)));
    Gov_spending_low_eta(t-start_time+1)     = eval(strcat('transition_traditional_fiscal_shocks.G',num2str(t)));  
end

%%%% financial shocks, Ramsey model with financial frictions
load transition_fiscal_shocks_eta_positive
for t = start_time:end_time+1  
    K_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.K',num2str(t)));
    B_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.B',num2str(t)));    
    Y_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.Y',num2str(t)));    
    B_Y_high_eta(t-start_time+1)   = eval(strcat('transition_fiscal_shocks.B_Y',num2str(t)));
    C_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.C',num2str(t)));
    L_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.L',num2str(t)));
    phi_high_eta(t-start_time+1)   = eval(strcat('transition_fiscal_shocks.phi',num2str(t)));
    I_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.I',num2str(t)));
    q_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.q',num2str(t)));
    tau_l_high_eta(t-start_time+1) = eval(strcat('transition_fiscal_shocks.tau_l',num2str(t)));
    tau_k_high_eta(t-start_time+1) = eval(strcat('transition_fiscal_shocks.tau_k',num2str(t)));
    R_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.R',num2str(t)));
    Gov_spending_high_eta(t-start_time+1)     = eval(strcat('transition_fiscal_shocks.G',num2str(t)));    
    q_w_high_eta(t-start_time+1)   = eval(strcat('transition_fiscal_shocks.q_w',num2str(t)));
    r_K_high_eta(t-start_time+1)   = eval(strcat('transition_fiscal_shocks.r_K',num2str(t)));
    rho_high_eta(t-start_time+1)   = eval(strcat('transition_fiscal_shocks.rho',num2str(t)));
end


%% Comparing without financial frictions and the other with financial frictions (always no shocks)
start_time = 0; 
end_time   = 20;
figure()
subplot(3,4,1)
h = plot(time_period, I_low_baseline ./ K_high_baseline(end) / par.DELTA * 100, time_period, I_high_baseline ./ K_high_baseline(end) / par.DELTA * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Investment')
ylabel('% of ss')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,2)
h = plot(time_period, C_low_baseline ./ C_high_baseline(end) * 100, time_period, C_high_baseline ./ C_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Consumption')
ylabel('% of ss')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,3)
h = plot(time_period, L_low_baseline  ./ L_high_baseline(end) * 100, time_period, L_high_baseline ./ L_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Hours')
ylabel('% of ss')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,4)
h = plot(time_period, Y_low_baseline ./ Y_high_baseline(end) * 100, time_period, Y_high_baseline ./ Y_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Output')
ylabel('% of ss')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,5)
h = plot(time_period, K_low_baseline ./ K_high_baseline(end) * 100, time_period, K_high_baseline ./ K_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('K')
ylabel('% of ss')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,6)
h = plot(time_period, B_low_baseline / B_high_baseline(end) * 100, time_period, B_high_baseline / B_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Debt')
ylabel('% of ss')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,7)
h = plot(time_period, phi_low_baseline * 100, time_period, phi_high_baseline * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$\phi$','interpreter','Latex')
ylabel('%')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,8)
h = plot(time_period, q_low_baseline, time_period, q_high_baseline);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$q$','interpreter','Latex')
ylabel('')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,9)
h = plot(time_period, tau_l_low_baseline * 100, time_period, tau_l_high_baseline * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$\tau^{\ell}$','interpreter','Latex')
ylabel('%')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,10)
temp=subplot(3,4,10);
temppos=get(temp,'position');
h = plot(time_period, tau_k_low_baseline * 100, time_period, tau_k_high_baseline * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$\tau^k$','interpreter','Latex')
ylabel('%')
set(gca, 'FontSize', 16)
xlim([start_time end_time])
ylim([0 50])
annotation('textarrow',[temppos(1)+.055*temppos(3) temppos(1)+.055*temppos(3)],[temppos(2)+.58*temppos(4) temppos(2)+.98*temppos(4)],'String',compose("%3.0f",tau_k_low_baseline(2)*100) + "%",'Fontsize',16)
annotation('textarrow',[temppos(1)+.065*temppos(3) temppos(1)+.065*temppos(3)],[temppos(2)+.26*temppos(4) temppos(2)+.98*temppos(4)],'String',compose("%3.0f",tau_k_high_baseline(2)*100) + "%",'Fontsize',16,'LineStyle','-.','color','r')
grid on

subplot(3,4,11)
h = plot(time_period, R_low_baseline * 100 - 100, time_period, R_high_baseline * 100 - 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Interest rate')
ylabel('%')
xlim([start_time end_time])
ylim([-50, 10])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,12)
h = plot(time_period, (Gov_spending_low_baseline./par.G - 1 )* 100, time_period, (Gov_spending_high_baseline /par.G -1 ) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$G$','interpreter','Latex')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on
legend('no financial frictions','financial frictions','Orientation','horizontal', 'Fontsize',14)



%% Comparing without financial frictions and the other with financial frictions (always with fiscal shocks)
start_time = 5; 
end_time = 25;
figure()
subplot(3,4,1)
h = plot(time_period, (I_low_eta - I_low_baseline) ./ K_high_baseline(end) / par.DELTA * 100, time_period, (I_high_eta - I_high_baseline) ./ K_high_baseline(end) / par.DELTA * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Investment')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,2)
h = plot(time_period, (C_low_eta - C_low_baseline) ./ C_high_baseline(end) * 100, time_period, (C_high_eta - C_high_baseline) ./ C_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Consumption')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,3)
h = plot(time_period, (L_low_eta - L_low_baseline ) ./ L_high_baseline(end) * 100, time_period, (L_high_eta - L_high_baseline) ./ L_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Hours')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on


subplot(3,4,4)
h = plot(time_period, (Y_low_eta - Y_low_baseline) ./ Y_high_baseline(end) * 100, time_period, (Y_high_eta - Y_high_baseline) ./ Y_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Output')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,5)
h = plot(time_period, (K_low_eta - K_low_baseline) ./ K_high_baseline(end) * 100, time_period, (K_high_eta - K_high_baseline) ./ K_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('K')
ylabel('% dev')
%ylim([-1,end_time])
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,6)
h = plot(time_period, (B_low_eta - B_low_baseline) / B_high_baseline(end) * 100, time_period, (B_high_eta - B_high_baseline) / B_high_baseline(end) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Debt')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,7)
h = plot(time_period, (phi_low_eta - phi_low_baseline) * 100, time_period, (phi_high_eta - phi_high_baseline) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$\phi$','interpreter','latex')
ylabel('diff in percent')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,8)
h = plot(time_period, (q_low_eta - q_low_baseline) * 100, time_period, (q_high_eta - q_high_baseline) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$q$','interpreter','latex')
ylabel('diff in percent')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,9)
h = plot(time_period, (tau_l_low_eta - tau_l_low_baseline) * 100, time_period, (tau_l_high_eta - tau_l_high_baseline) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$\tau^{\ell}$','interpreter','latex')
ylabel('diff in percent')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,10)
h = plot(time_period, (tau_k_low_eta - tau_k_low_baseline) * 100, time_period, (tau_k_high_eta  - tau_k_high_baseline) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$\tau^k$','interpreter','latex')
ylabel('diff in percent')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,11)
h = plot(time_period, (R_low_eta - R_low_baseline) * 100, time_period, (R_high_eta - R_high_baseline) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('Interest rate')
ylabel('diff in percent')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on

subplot(3,4,12)
h = plot(time_period, (Gov_spending_low_eta./par.G - 1 )* 100, time_period, (Gov_spending_high_eta /par.G -1 ) * 100);
set(h(1),'LineStyle','-','Linewidth',2,'color','k')
set(h(2),'LineStyle','-.','Linewidth',2,'color','r')
title('$G$','interpreter','Latex')
ylabel('% dev')
xlim([start_time end_time])
set(gca, 'FontSize', 16)
grid on
legend('no financial frictions','financial frictions',...
    'location','best', 'Orientation','horizontal', 'Fontsize',14)



%% Check residuals of Euler equations and goods-market clearing conditions for fiscal shocks with financial frictions (optional)
display('Done...Well, I am going to do a final check of the equilibrium objects...')
display('Just to be sure we get things right...')
display(' ')

end_period      = length(tau_k_high_eta) - 1; % starting from period 0
MRS             = par.BETA * C_high_eta(1:end-1) ./ C_high_eta(2:end);
tau_k_plus_1    = tau_k_high_eta(2:end);
r_K_plus_1      = r_K_high_eta(2:end);
phi_plus_1      = phi_high_eta(2:end);
q_plus_1        = q_high_eta(2:end);
rho_plus_1      = rho_high_eta(2:end);
q_w             = q_w_high_eta(1:end-1);
q_w_plus_1      = q_w_high_eta(2:end);
time_           = 0:(end_period - 1);

Consumption     = C_high_eta(2:end);
Government      = Gov_spending_high_eta(2:end);
Capital         = K_high_eta(2:end);
Labor           = L_high_eta(2:end);
Capital_minus_1 = K_high_eta(1:end-1);
liquidity       = phi_high_eta(2:end);
int_cost        = par.eta_0 * liquidity.^(par.eta_1);
time_resources  = 1:end_period;

rid_resources = Consumption + Government + Capital .* (1 + liquidity .* int_cost) - ...
    (...
    int_cost .* liquidity .* (1 - par.CHI) .* (1 - par.DELTA) .* Capital_minus_1 + par.A * Capital_minus_1.^(par.ALPHA) .* Labor.^(1 - par.ALPHA)...
    + (1 - par.DELTA) .* Capital_minus_1...
    );

rid_EE_capital             = q_w  - MRS .*... 
    (...
    r_K_plus_1 .* (1 - tau_k_plus_1) + tau_k_plus_1 .* q_w_plus_1 * par.DELTA + (1 - par.DELTA) * q_w_plus_1 + par.CHI * (1 - par.DELTA) * phi_plus_1 .* (q_plus_1 .* (1 + rho_plus_1) - q_w_plus_1)...
    );
display('time, q_w, MRS, capital tax, r_K, phi, q, rho, Euler residual')

matrix_EE_capital  = [time_', q_w_high_eta(1:end-1), MRS, tau_k_plus_1, r_K_plus_1, phi_plus_1, q_plus_1, rho_plus_1, rid_EE_capital]

display('time, C, G, K_minus_1, L, phi, eta, Goods mkt residual')
matrix_resoucres   = [time_resources', Consumption, Government, Capital_minus_1, Labor, liquidity, int_cost, rid_resources]
    










